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Confinement breaks the translational symmetry of materials, making all thermodynamic and ki¬ 
netic quantities functions of position. Such symmetry breaking can be used to obtain conhgurations 
that are not otherwise accessible in the bulk. Here, we use computer simulations to explore the effect 
of substrate-liquid interactions on thermodynamic and kinetic anisotropies induced by a solid sub¬ 
strate. We consider n-octane nano hlms that are in contact with substrates with varying degrees of 
attraction, parameterized by an interaction parameter es- Complete freezing of octane nano hlms 
is observed at low temperatures, irrespective of e^, while at intermediate temperatures, a frozen 
monolayer emerges at solid-liquid and vapor-liquid interfaces. By carefully inspecting the prohles of 
translational and orientational relaxation times, we conhrm that the translational and orientational 
degrees of freedom are decoupled at these frozen monolayers. At sufficiently high temperatures, 
however, free interfaces and solid-liquid interfaces close to loose (low-e^) substrates undergo ’pre- 
freezing‘, characterized by mild peaks in several thermodynamic quantities. Two distinct dynamic 
regimes are observed at solid-liquid interfaces. The dynamics is accelerated in the vicinity of loose 
substrates, while sticky (high-e^) substrates decelerate dynamics, sometimes by as much as two 
orders of magnitude. These two distinct dynamical regimes have been previously reported by us 
[J. Chem. Phys. 141 : 024506 (2014)] for a model atomic glass-forming liquid. We also conhrm 
the existence of two correlations- proposed in the above-mentioned work- in solid-liquid subsur¬ 
face regions of octane thin hlms, i.e., a correlation between atomic density and normal stress, and 
between atomic translational relaxation time and lateral stress. Finally, we inspect the ability of 
different regions of an octane him to explore the potential energy landscape by performing inherent 
structure calculations, and observe no noticeable difference between the free surface and the bulk in 
efficiently exploring the potential energy landscape. This is unlike the hlms of model atomic glass 
formers that tend to sample their respective landscape more efficiently at free surfaces. We discuss 
the implications of this hnding to the ability of octane- and other n-alkanes- to form ultrastable 
glasses. 


I. INTRODUCTION 

Conhnement alters the thermodynamic and kinetic 
properties of matter. Such changes can partly arise from 
quantum effects at the nanoscale [1-4]. But they can 
also emerge in purely classical systems, simply due to 
the presence of conhnement-induced interfacial regions. 
On a fundamental level, conhnement breaks the transla¬ 
tional isotropy of a bulk material, making all its phys¬ 
ical properties functions of position [5]. The extent of 
such anisotropy is, however, variable, and can depend on 
both the thermodynamic conditions (temperature, den¬ 
sity, composition) as well as the nature of the interface(s), 
such as boundary conditions (free interface vs. substrate), 
interactions (attractive, repulsive, inert) and dimension¬ 
ality (pores, channels, slabs). Regardless, the behavior of 
conhned matter can deviate signihcantly from the bulk in 
many different ways. Examples include the emergence of 
new phases [6-8], changes in thermodynamic properties 
such as melting and glass transition temperatures [9-13], 
and in kinetic properties such as nucleation rates [14- 
16], viscosities [17-20], diffusivities [21-23] and elastic 
constants [18, 24, 25]. 

The translational anisotropy that is induced as a result 
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of conhnement can be utilized for modulating the struc¬ 
tural and functional properties of materials. A widely 
known example is heterogeneous nucleation [26] in which 
a substrate enhances a particular hrst-order phase tran¬ 
sition in its vicinity by decreasing the associated free en¬ 
ergy barriers. On a microscopic level, this is mediated 
by the formation of a subsurface region that is more con¬ 
ducive to nucleation than the homogeneous bulk mate¬ 
rial. The nucleating potency of a substrate is thus related 
to the microstructure of the material that is in its proxim¬ 
ity, i.e., the anisotropy that it induces therein [27]. The 
relationship between the anisotropy of the subsurface re¬ 
gion and the facilitated nucleation can be rigorously ex¬ 
plained for phase transitions such as hydrophobic evapo¬ 
ration [14, 28], but is far more difficult to discern for more 
complex phase transitions, such as crystallization [16]. 

A more recent example of the interesting behavior 
arising in anisotropic systems is provided by ultrastable 
glasses obtained by depositing the vapor of a glass¬ 
forming liquid onto a cold substrate. By tuning the sub¬ 
strate temperature, it is possible to form glasses that are 
far more stable than the ordinary glasses obtained by 
rapidly quenching the liquid [29] . In addition to their su¬ 
perior thermal and mechanical stability, such ultrastable 
glasses can be structurally anisotropic as evident from 
refractive index measurements and X-ray scattering [30- 
32] . Possible structural differences between ordinary and 
ultrastable glasses are manifest in the different scaling of 
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their heat capacities in the limit of T -> 0 K [33]. It 
has even been suggested that a first-order liquid-liquid 
transition might exist between the ordinary and ultra¬ 
stable amorphous states [34]. Vapor deposition has been 
successfully used for making organic [29, 35-41], poly¬ 
meric [42] and metallic [43] ultrastable glasses. There is, 
however, a significant gap in understanding why these 
vapor-deposited glasses are ultrastable. It has been ar¬ 
gued that the enhanced mobility of molecules at the 
vapor-liquid interface enables them to explore the po¬ 
tential energy landscape more efficiently, thereby giving 
rise to structures that reside far deeper in the poten¬ 
tial energy landscape [44]. Such enhanced mobility has 
been observed in experimental [20, 45-47] and compu¬ 
tational [5, 48-50] studies of thin films. However, no 
causal relationship has been unambiguously established 
between the existence of such accelerated regions and the 
formation of ultrastable glasses 

A systematic way of identifying a possible link between 
enhanced surface mobility and increased thermodynamic 
and kinetic stability is to probe thermodynamic and ki¬ 
netic anisotropies in films of different compositions, and 
to assess their sensitivity to changes in thermodynamic 
conditions. It is of particular interest to study the role 
of substrates in inducing such anisotropies, and in mod¬ 
ifying the properties of the resulting glasses. This is 
a question that has not been extensively studied, espe¬ 
cially in regard to the vapor-deposited ultrastable glasses. 
Two particular parameters that are relevant in this quest 
are the temperature and the interaction parameter, with 
the latter quantifying the strength of (attractive) inter¬ 
actions between the substrate and the liquid. In our ear¬ 
lier publication [5], we performed a detailed analysis of 
such anisotropies in thin films of a model atomic glass¬ 
forming liquid [51], and we discovered two distinct dy¬ 
namical regimes: accelerated dynamics near loosely at¬ 
tractive substrates, and decelerated dynamics near sticky 
substrates. We also established correlations between os¬ 
cillations in density and in normal stress, and between 
oscillations in relaxation time and in lateral stress. 

In this work, we revisit the main findings of Ref. [5], 
but now in the context of molecular thin films of n- 
octane. We choose n-octane as a prototypical chain 
molecule. Alkane thin films are particularly interesting 
as they can undergo a process known as surface freez¬ 
ing [52, 53] in which a frozen layer emerges at the vapor- 
liquid interface at temperatures exceeding the equilib¬ 
rium melting temperature. Also, alkanes constitute one 
of the most important components of crude oil, and their 
presence can result in interesting phase separation phe¬ 
nomena that are typically governed by confinement [54- 
57]. It is therefore worthwhile to inspect the microstruc¬ 
ture of interfacial regions of alkane films with an eye to¬ 
wards understanding the mechanism(s) and predicting 
the kinetics of such phase transition processes. 

This paper is organized as follows. Section IIA pro¬ 
vides details on the computational setup as well as the 
force-fields employed. Technical details of the Simula- 





FIG. 1: (a) Representation of an n-octane molecule with the 
NERD force field. Green and light bone sites correspond to 
the united-atom GH 3 and GH 2 interaction sites respectively, 
(b) A liquid film of n-octane molecules in the vicinity of an 
attractive (steel blue) substrate. The light pink substrate is 
repulsive. 


tions performed are given in Section IIB. The proce¬ 
dures utilized for computing spatial profiles of thermo¬ 
dynamic and kinetic properties are given in Section IIC. 
Section III A discusses the qualitative behavior of the 
octane films, particularly, in regard to surface freezing. 
Anisotropies in thermodynamic and kinetic properties 
are discussed in Sections IIIB and IIIC, respectively. Fi¬ 
nally, Section IV is reserved for concluding remarks. 


II. METHODS 
A. System Description 

A schematic representation of the n-octane films con¬ 
sidered in this work is depicted in Fig. lb. Octane 
molecules are represented with the NERD force-field [58] 
in which each octane molecule is comprised of eight 
united-atom interaction sites: two for CH 3 (green in 
Fig. la) and six for CH 2 (light blue in Fig. la). Our 
choice of a united-atom potential such as NERD is due 
to its simplicity, as we only consider n-alkanes as pro¬ 
totypical chain molecules. Indeed, the increased quanti¬ 
tative accuracy obtained from using more realistic, but 
computationally costly, all-atom force-fields is unlikely 
to impact the main findings of this work, and can only 
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shift the n values or temperatures for which the reported 
phenomena are observed. In the NERD potential, the 
non-bonded interactions between these sites are modeled 
through the Lennard-Jones potential [59]: 


VLj{r) = 4ey 



The bonded interactions include bond-stretching, bond¬ 
bending and torsional terms: 


Stretching ) 
Sending (^) 


Sorsional (0) 


ksir-rof (2a) 

h{e-0of (2b) 

+ COS0) + ^^ 2(1 + cos 20) 


+ cos 30) 


(2c) 


All model parameters are given in Table. I. Lennard- 
Jones interactions are shifted and truncated at Vc = 
1.38 nm. 

In addition to the n-octane film, there are two sub¬ 
strates in the system: the attractive substrate {C atoms, 
steel blue in Fig. lb) and the repulsive substrate {D 
atoms, light pink in Fig. lb). The constituent atoms 
of both substrates are arranged into a face-centered cu¬ 
bic (fee) lattice. They interact with themselves and the 
united-atom LJ sites via the Lennard-Jones potential 
with cac = ^AD = ^s^AAi ^BC = ^BD = ^s^BB, and 
O'AC = obc = oad = obd = 0.4 nm. Here, the in¬ 
teraction parameter, e^, is used for tuning the strength 
of attractive (repulsive) interactions between the sub¬ 
strates and octane molecules. All substrate-octane in¬ 
teractions are truncated and shifted at Vc = 1.2 nm and 
Vc =0.44 nm for the attractive and repulsive substrates, 
respectively. The inclusion of a second repulsive sub¬ 
strate is to assure that the octane molecules evaporating 
from the film would never redeposit onto the opposite 
side of the substrate as a result of periodic boundary 
conditions. Simulations are carried out for three distinct 
values of es = 0.5,1.0 and 3.0, with the temperature 
range 200 K < T < 290 K considered at each es value. 
As will be shown in Section III, this temperature range 
covers the states that are both above and below freezing. 


B. Simulation Details 

All Molecular Dynamics (MD) simulations are per¬ 
formed using LAMMPS [60] in the isochoric (NVT) en¬ 
semble. We integrate Newton’s equations of motion us¬ 
ing the velocity Verlet algorithm [61] with a time step of 
2 fs, and we control temperature using the Nose-Hoover 
thermostat [62, 63] with a time constant of r = 0.2 ps. 
All simulations are carried out in cuboidal boxes that are 
periodic in all dimensions. The simulation box is always 
longer along the z direction in order to assure the lack 
of correlation between the liquid film and its periodic 
images. 


TABLE I: Parameters of the NERD force-field for n-octane 
molecules. A and B correspond to the united-atom CH 3 and 
CH 2 sites, respectively. 

Parameter Value 


eAA 
cr A A 
esB 
ctbb 

eAB = \/eAAeBB 

CTAB = (ctaa + crsn)/2 

ks 

ro 

kb 

Oo 

Ki 

K 2 

Ks 


0.2066 kcal/mol 
0.391 nm 
0.09106 kcal/mol 
0.393 nm 
0.1372 kcal/mol 
0.392 nm 

191.85 kcal/mol-A^ 
0.154 nm 

124.26 kcal/mol-radian' 
114° 

1.4117 kcal/mol 
—0.2711 kcal/mol 
3.1465 kcal/mol 


In all simulations, the starting configuration is com¬ 
prised of an n-octane film in which all molecules are ar¬ 
ranged into a simple cubic lattice. This configuration is 
initially heated at T = 300 K for 200 ps in order to melt 
the crystal. The melted film is then gradually quenched 
to Tj at a cooling rate of 2.5 x 10^^ K/s. The arising 
configuration is equilibrated at Tj for 4 ns. This equili¬ 
bration time is far larger than the structural relaxation 
time in the bulk for the temperatures considered in this 
work. Throughout the entire process, a separate ther¬ 
mostat is applied to the atoms in the substrate, always 
maintained at a temperature Tf. After this initial equi¬ 
libration stage, the production runs are carried out for 
80 ns. 


C. Spatial Profiles 

1. Thermodynamic Properties 

The bulk of the methodology that is used for com¬ 
puting the spatial profiles of thermodynamic and kinetic 
properties is discussed in our earlier publication [5]. For 
thermodynamic quantities that are time-invariant, spa¬ 
tial profiles can be rigorously determined via simple time 
averaging of the corresponding property in thin cuboidal 
slices of the simulation box. Here, profiles of potential 
energy, atomic and molecular density, lateral and normal 
stress and lateral radial distribution function are com¬ 
puted in slabs that are 0.025 nm thick, and the binning 
of all molecular properties, including molecular density, 
and radial distribution function, is performed based on 
the centers of mass of the molecules. We also compute 
inherent structures using the FIRE algorithm [64], and 
the contribution to the average inherent structure poten¬ 
tial energy of a slice is from those atoms that were orig- 
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inally located in that slice prior to energy minimization. 
We also compute the orientational distribution function 
(ODF), f{0^z), of n-octane molecules defined as follows. 
First, % the gyration tensor of each n-octane molecule is 
computed from: 


= 






( 3 ) 




Here rriij and Vij are the mass and the position of the 
jf'th united atom site of the ith molecule and Vi^cM is the 
center of mass of the ith molecule. Subsequently, v^, the 
longest principal axis of the ith molecule is determined, 
which is the unit eigenvector corresponding to the largest 
eigenvalue of Due to the inversion symmetry of 
both will be equally valid choices. In order to avoid 
any ambiguity, we choose the for which > 0. 

Here, is a unit vector perpendicular to the substrate, 
and pointing towards the liquid. The orientational dis¬ 
tribution function is then defined as: 




(Eili'^dv* -nsl - cose) 5{zi^cM - z)'j 

{'EIiS{z,,cm-z)) 


( 4 ) 


Note that jQ^^f{0^z)sm0d0 = 1. f{0,z) is utilized for 
computing further orientational order parameters that 
are introduced and discussed in Section HI A. 


2. Kinetic Properties 

The time averaging procedure that is used for comput¬ 
ing profiles of thermodynamic properties cannot be uti¬ 
lized for kinetic properties such as relaxation times as the 
latter are obtained from ensemble averages of autocorre¬ 
lation functions. This difficulty stems from the ambiguity 
of defining autocorrelation functions in open systems. In 
our earlier publication [5] , we present a thorough discus¬ 
sion of different heuristics that can be used for defining 
autocorrelation functions in confined systems. Here, we 
take the same convention as used in that work, and de¬ 
fine atomic and molecular translational and rotational 
relaxation times as follows. For translation relaxation, 
we compute the ^-dependent self intermediate scattering 
function: 


Fs,x{q, z,t) 




( 5 ) 


the beginning and at the end of the time window, and 
xVx{z^t) = Zi{ff)^ z)) is the average num¬ 

ber of entities that contribute to the sum for a particular 
slice in Eq. (5). Like thermodynamic properties, contri¬ 
butions to molecular autocorrelation functions are based 
on the centers of mass of n-octane molecules. For oc¬ 
tane films, we use a g = 16 nm“^, which corresponds 
to the first peak of 5'(g), the structure factor, in the 
bulk n-octane liquid computed in an NpT simulation at 
T = 300 K and p = 0 bar. We observe no significant 
change in the locus of the maximum of S{q) in bulk oc¬ 
tane with temperature. Relaxation times are determined 
from Fs,x{q, z^t = r) = 0.2. The atomic and molecular 
translational relaxation times are denoted by rtr,a and 
^tr,m, respectively. 

In order to quantify molecular rotational relaxation 
time profiles, we compute the following ^-dependent ori¬ 
entational auto-correlation function: 


/Nm \ 


(6) 


with P 2 (cos = (3/2) cos^ ^ —(1/2) the second Legendre 
polynomial. Note that P 2 (') is the natural choice consid¬ 
ering the inversion symmetry of the gyration tensor used 
for determining v^. Analogously, the relaxation of orien¬ 
tations of individual bonds is characterized using hi){z^ t), 
the bond autocorrelation function, defined as follows: 




hb{z,t) = 


Ji{z,t) 


^Pi[hi{t) • hi{0)]^{zi{t),zi{0),z) 


\ i=l 


( 7 ) 


Here is the unit vector in the direction of bond i that 
connects two united-atom sites on an n-octane molecule, 
and Zi is the 2 : coordinate of the center of the bond. Since 
a bond is essentially a directed entity with no inversion 
symmetry, we use Pi(cos6>) = cosO, the first Legendre 
polynomial, to quantify its rotational relaxation, as us¬ 
ing P 2 (') will lead to loss of relevant orientational infor¬ 
mation. Similar to translation relaxation times, molecu¬ 
lar and bond rotational relaxation times are determined 
from hm{z,t = = 0.2 and hi,{z,t = Trot.b) = 0.2, 

respectively. 

Since it is more difficult to obtain suitable statistics 
for the 2 :-dependent autocorrelation functions introduced 
above, we use slices that are 0.1 nm thick, four times 
thicker than the slices used for computing ^-dependent 
time averages of thermodynamic properties 


III. RESULTS AND DISCUSSION 


Here X = a, m, with a and m corresponding to atomic 
and molecular self-intermediate scattering functions, re¬ 
spectively. Ar (t) is the lateral displacement of entity i 
over time t. "Dizi, Z 2 ^ z) = S{zi — z)S(z 2 — z) is an in¬ 
dicator that assures that entity i is present at z both in 


A. Qualitative Behavior of Octane Films 

As mentioned in Section I, sufficiently long n-alkanes 
can undergo a process known as surface freezing [52, 53]. 
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FIG. 2: (a) Complete freezing of an n-octane film at 220 K, (b-c) Surface freezing of two n-octane films at 240 K in the vicinity 
of substrates with (b) es' = 0.5 and (c) es — 3.0. Lateral radial distribution functions reveal long-range translation order across 
the film at 220 K, while the 240-K films are only translationally ordered at the solid-liquid and vapor-liquid interfaces. 


To be more precise, alkane films behave differently at 
temperatures above and below T^, their effective melt¬ 
ing temperature. For films that are sufficiently thick 
Tm ~ the equilibrium bulk melting temperature. At 
T < T^, the entire film freezes into a lamellar crystal. At 
Tm < T < Tg, however, only a frozen monolayer emerges 
at the vapor-liquid interface. Here stands for the sur¬ 
face freezing temperature. Experimentally, surface freez¬ 
ing occurs for 15 < n < 50 [65]. Yet, we observe sur¬ 
face freezing in n-octane films simulated using the NERD 
force-field, even though n = 8 for n-octane. Eor the 
films considered in this work, complete freezing occurs 
for T < 230 K irrespective of the es value (Eig. 2a). The 
existence of long-range lateral order in these ’cold‘ films is 
clearly visible in the lateral radial distribution functions 
depicted in Eig. 2a. Eor 235 K < T < 250 K, however, 
only a frozen monolayer emerges at the vapor-liquid and 
solid-liquid interfaces, with the center of the film remain¬ 
ing amorphous (Eigs. 2b-c). Here, we define freezing in 
a purely structural sense, i.e., based on the presence of 
long-range lateral translational order. As will be fur¬ 
ther elucidated in Section HIC, such frozen regions also 
correspond to dynamically decelerated regions in which 
translational and rotational degrees of freedom are de¬ 
coupled. 

The crossover temperature of Tm = 232.5 ±2.5 K, 
is reasonably close to 216 K, the experimental melting 
temperature of n-octane [66]. (No computational esti¬ 
mate of Tm is available for the NERD force-field.) Our 
findings are also qualitatively consistent with earlier com¬ 
putational studies of n-octane [67, 68], and other closely- 
related n-alkanes [69] using different force-fields, all re¬ 
porting surface freezing. 

In order to examine the microstructure of frozen mono¬ 


layers, we compute the ODEs defined in Eq. (4). Eig. 3 
depicts the ODEs for films at T = 240 K. The frozen 
monolayers emerging at vapor-liquid interfaces are char¬ 
acterized by a single peak in /(6 >,z), corresponding to 
a perpendicular arrangement of octane molecules. Such 
an arrangement is easily visible in the films depicted in 
Eig. 2a-c. In the vicinity of the substrate, however, the 
microstructure of the frozen monolayer depends on e^. 
Eor loose substrates, i.e. es = 0.5 (Eig. 3a) and = 1 
(Eig. 3b), the microstructure closely resembles that of 
the free interface, as evident in the single peak of /(6>, z) 
at 6> = 0 and z ^ 1 nm. In contrast, a sticky sub¬ 
strate induces a different type of ordering, with the aris¬ 
ing ODE having two distinct strong peaks at ^ = 0: one 
at z = 0.65 nm, and one at 2 ; = 0.9 nm. The loci of these 
peaks are identical to the peaks in the molecular density 
profile depicted in the rightmost panel of Eig. 6. The 
existence of two peaks in /(^, z) is a consequence of the 
corrugated surface of the 001 facet of fee substrate. In 
other words, the molecules at the frozen monolayer can 
be present both in the valleys and peaks of the rough 001 
surface. This changes their z value even though they all 
have the same orientation with ^ = 0. 

At higher temperatures, no frozen monolayer emerges 
at vapor-liquid interfaces (Eig. 4b). However, n-octane 
molecules tend to have a mild preponderance to align 
along the 2 : axis, as evident in the ODE depicted in 
Eig. 4a. Such a propensity can be qualitatively described 
as ’pre-freezing‘, a phenomenon previously observed in 
earlier computational studies of alkane films [68]. The 
notion of pre-freezing here is different from and should 
not be confused with the pre-freezing that is occasionally 
used for describing the interfacial freezing that precedes 
bulk freezing at T > Tm [70]. We can quantify the ex¬ 
tent for such pre-freezing with the following two order 
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FIG. 3: Orientational distribution functions for a film at 
T = 240 K and (a) es = 0.5, (b) = 1.0, (c) = 3.0. 


parameters. The first one is ^ = max; 2 > 6 nm /(^, the 
maximum of the free-surface peak of ODF. Note that 
^ = 1 for a fully isotropic film. Any deviation from unity 
will henceforth correspond to broken rotational symme¬ 
try. As depicted in Fig. 4c, ^ is always significantly larger 
than unity, even for the films at T = 290 K. Also the peak 
always occurs at = argmax^> g nm /(^5 = 0 irrespec¬ 

tive of T. As expected, pre-freezing becomes stronger at 
lower temperatures. 

The second order parameter that is adopted from the 
liquid crystal literature is called the nematic order pa¬ 
rameter (OP) and is essentially the second moment of 
the ODF [71]: 


Szz{z) = ^ / [3cos^ 6> — l]/(6 >, z) sin6>d6> (8) 

2 Jo 

For a fully isotropic fluid, Szz = 0, while the values of 
+1 and —0.5 correspond to perfect alignment along and 
perpendicular to the 2 : axis, respectively. Fig. 5 depicts 
Szz{z) profiles for different films. Note that the nematic 
OP is very close to unity for surface-frozen monolay¬ 
ers at low temperatures. Surface pre-freezing of high- 
temperature films is manifest in positive peaks of the ne- 




Temperature [K] 


FIG. 4: (a) Orientational distribution function at T = 270 K 
and es — 0.5, (b) A characteristic snapshot of the correspond¬ 
ing octane film. No freezing is observed, (c) Free surface 
alignment propensities at high temperatures. 


matic OP corresponding to weak alignment along the z 
direction. Similar to ^ that decreases with T, the heights 
of these peaks diminish as T increases. 

Similar to the free interface, a mildly pre-frozen mono- 
layer emerges in the vicinity of loose substrates {es = 
0.5,1.0), with characteristic mild peaks in the nematic 
OP profiles (Fig. 5). For sticky substrates, however, a 
frozen monolayer emerges at all temperatures, with its 
ODF resembling the one depicted in Fig. 3c. The nematic 
order parameter also demonstrates two strong peaks close 
to unity, consistent with perfect alignment along the 2 : di¬ 
rection spatially modulated by the corrugated surface of 
the substrate. 


B. Thermodynamic Properties 

Density: Profiles of atomic and molecular densities are 
depicted in Fig. 6. Both profiles are peaked at the frozen 
monolayers that emerge at the vapor-liquid and solid- 
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FIG. 5: Profiles of nematic order parameter for different temperatures and values. 
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FIG. 6: Profiles of atomic and molecular density for different temperatures and es values. 


liquid interfaces. In the case of free interfaces and loosely 
attractive substrates, a single large peak is observed. For 
sticky substrates, however, the peak splits into two as a 
result of spatial modulation induced by the corrugated 
substrate. The loci of these peaks correspond to the 
maxima of ODF (Fig. 4) and nematic OP (Fig. 5). Such 
large peaks are always followed by distinct valleys that 
correspond to the crystal/liquid interface. As evident in 
Fig. 5, the nematic OP is negative at these valleys as 
the octane molecules have a (weak) propensity to align 
parallel to the substrate. 

A very interesting feature of octane films is the strati¬ 
fication of octane molecules at free interfaces even in the 
absence of surface freezing, i.e., at T > 250 K. This be¬ 
havior is contrary to what is observed in atomic thin films 
in which density drops monotonically across a free inter¬ 
face [5, 48]. A single mild peak emerges in the atomic and 
molecular density profile, with its amplitude closely cor¬ 
relating with the pre-freezing order parameters defined in 
Section III A. A linear correlation between ^ (Fig. 4c) and 
Pmax,free 5 the amplitude of the molecular density peak, is 
depicted in Fig. 7. Henceforth, the nonmonotonic behav¬ 
ior of density is a result of pre-freezing that creates local 
high-density regions in the interfacial region. Density 
oscillations at a free interface have been previously ob¬ 
served in experimental and computational studies of not 



FIG. 7: Linear correlation between ^ and pmax,free. The dark 
line is a linear fit to the individual data points. 


only alkanes [68, 72] but also alkali metals [73, 74] and 
ionic liquids [75]. The emergence of density oscillations 
has been attributed to a wide range of factors, includ¬ 
ing large separation between the critical and triple-point 
temperatures [76] and building block anisotropy [77]. Un¬ 
like the atomic thin films, the oscillatory density profiles 
in octane films cannot be fitted to the commonly used 
hyperbolic tangent functional form [78], henceforth mak¬ 
ing the determination of the width of the free interface 





































based on density nontrivial. 

Close to the substrates, both atomic and molecular 
density profiles are oscillatory. This is consistent with 
the traditional picture of confinement in which a wall 
induces structure in the surrounding liquid [79-81]. The 
stratification of the octane liquid is, however, an interplay 
between the structuring induced by the substrate and the 
natural propensity of the interfacial region to pre-freeze. 
Density oscillations in the vicinity of a loosely attractive 
substrate {es = 0.5) are very similar to the profiles at 
the free interface, with the substrate inducing very little 
structure in the liquid. This is very similar to theoretical 
predictions obtained from density functional theory for 
associating fluids near a hard wall [81]. As es increases, 
however, stratification becomes more pronounced. For 
instance, multiple peaks in density emerge at es = 1.0, 
even though the solid-liquid interface has the same ori¬ 
entational fingerprints as the vapor-liquid interface. 
Potential Energy and Internal Structure: Profiles 
of (t/), the average potential energy, are depicted in 
Fig. 8. (U) exhibits a maximum at the center of a frozen 
monolayer. This apparent energetic penalty can be at¬ 
tributed to the lamellar microstructure of a monolayer in 
which the united-atom CH 2 sites cluster in the middle. 
Since |ecH 2 ,CH 2 l < kcH 2 ,CH 3 l < |ecH 3 ,CH 3 |, such CH 2 - 
rich regions will correspond to higher potential energies. 
This energetic penalty is, however, offset by the CHs-rich 
regions at two ends of a monolayer, corresponding to min¬ 
ima in (U). It is also noteworthy that z*, the locus of the 
maximum of (I/), matches the corresponding for den¬ 
sity, f{0,z) and Szz{z)^ further confirming that all these 
oscillations are manifestations of the same phenomenon, 
i.e., interfacial freezing. Like density, (U) profiles are os¬ 
cillatory close to the substrate, even in the absence of 
interfacial freezing, with the amplitudes and number of 
oscillations increasing with 1/T and e^. However, we do 
not observe any nonmonotonicity in the pre-frozen free 
interfaces, suggesting that the propensity of a surface to 
pre-freeze has a purely entropic origin. 

Although potential energy profiles provide valuable in¬ 
formation about the relative energetic stability of differ¬ 
ent regions of a confined material, they convey very little 
about the ability of a material to explore its confinement- 
induced potential energy landscape. A systematic way of 
characterizing the latter is to obtain an ensemble of inher¬ 
ent structures, and to compute the amount of decrease in 
potential energy for different regions of the material [48] . 
Here, 2 : binning is performed based on the initial (pre¬ 
minimization) positions of the individual atoms. Spatial 
profiles of this quantity, which we denoted as dive pro¬ 
files in our earlier publication [5], are depicted in Fig. 9. 
For all the films, whether they go through surface freez¬ 
ing or pre-freezing, no noticeable overall difference is de¬ 
tected between the landscape depth in the bulk and at 
the free interface. The ability of a free interface to modify 
the potential energy landscape accessible to a liquid is, 
therefore, far more limited in the case of n-octane. This 
is unlike the behavior observed in the atomic thin films 


considered in Refs. [5, 48] in which a vapor-liquid inter¬ 
face increases the depth of the potential energy landscape 
accessible to the liquid at its vicinity. Such a heightened 
access to the minima of the potential energy landscape 
appears to be key in the ability of a material to form 
ultrastable glasses upon vapor deposition [82]. By this 
token, the ability of n-alkanes to form ultrastable glasses 
will be limited considering the lack of elevated access to 
the potential energy landscape of the liquid at the free 
interface. 

Stress: Of all the thermodynamic quantities considered 
in this work, none oscillates more strongly than normal 
and lateral stress (Fig. 10). Such oscillations are particu¬ 
larly strong in frozen monolayers at both interfaces, and 
extend well beyond the frozen region. Indeed, sharp in¬ 
terfaces are maintained as a result of surface tension that 
stems directly from anisotropies of the stress tensor [83] . 
In the case of frozen monolayers, in particular, two such 
sharp interfaces are present: one between the substrate 
(or vapor) and the frozen monolayer, and the other be¬ 
tween the monolayer and the liquid. Consequently, oscil¬ 
lations in stress occur across a relatively wider region of 
the film. 

It can be argued that stress anisotropy is among the 
most systematic ways of defining what constitutes an in¬ 
terfacial region. Indeed, a region with anisotropic stress 
responds differently to mechanical stress, and will thus 
have mechanical properties distinct from the bulk. In this 
work, we define the bulk region as the largest connected 
region in which the difference between normal and lateral 
stress does not exceed 20 bar. Choosing this cutoff, which 
is always less than 10 per cent of the largest difference 
between lateral and normal stress at T =290 K, is neces¬ 
sary considering the difficulty of converging virial-based 
estimates of stress in molecular simulations. Applying 
this criterion yields an interfacial width that increases 
with temperature, from ^2.30 nm at 260 K to «2.45 nm 
at 290 K. The interfacial width, as determined from the 
anisotropy in stress tensor, is very difficult to measure 
in experiments. Nevertheless, the temperature scaling of 
the dynamical length scales in thin films of polystyrene 
glasses reveals a similar dependence of temperature, with 
the width of the highly mobile region increasing with 
temperature [84] . It can indeed be argued that a correla¬ 
tion must exist between the thermodynamic length scale 
(obtained from stress tensor) and the dynamical length 
scales (defined as per mobility). 

In our earlier work [5], we established a strong corre¬ 
lation between the substrate-induced profiles of normal 
density and normal stress in amorphous regions of atomic 
thin films. By inspecting the density and normal stress 
profiles, we explore the existence of such a correlation 
in octane films. Here, the quantity of relevance is the 
atomic- and not the molecular- density, as the positions 
of individual molecules are not sufficient for determin¬ 
ing the stress tensor. We observe a correlation between 
atomic density and normal stress in amorphous regions of 
octane films. For loosely attractive substrates (Fig. 11a), 
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FIG. 8: Profiles of potential energy for different temperatures and es values. 
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FIG. 9: Profiles of dive energy for different temperatures and es values. 


the correlation is weaker to the extent that the first max¬ 
imum in normal stress only corresponds to a shoulder in 
atomic density. In other words, the solid-liquid inter¬ 
facial region is very similar to a free interface, as the 
substrate exerts minimal effect in ordering and stratify¬ 
ing the liquid in its vicinity. In the case of stickier sub¬ 
strates (Fig. 11b), however, the correlation is much more 
pronounced, and the peaks and valleys of atomic den¬ 
sity and normal stress follow one another more closely. 
It is, of course, necessary to emphasize the absence of 
such correlations in parts of the film that have long-range 
translational order, as crystals can respond to stress in 
nontrivial ways. 


C. Kinetic Properties 

Figs. 12 and 13 depict profiles of translational and ro¬ 
tational relaxation times. In surface-frozen monolayers, a 
sharp increase in relaxation times is observed, with both 
the atomic and molecular relaxations times increasing by 
at least one order of magnitude with respect to their bulk 
values. We observe an even larger increase in molecular 
and bond orientational relaxation times. We are unable 
to quantify the extent of such an increase as the orien¬ 
tational autocorrelation functions given in Eqs. (6) and 
(7) never relax to zero in the timescale of our simula¬ 
tions. This decoupling of translational and orientational 
degrees of freedom is a direct consequence of layering at 
the surface, and is common in phases with long-range 


orientational order such as liquid crystals [85]. In such 
systems, the relaxation of translational degrees of free¬ 
dom proceeds through a slipping mechanism that does 
not require orientational relaxation. 

In Section IIIB, we discuss non-monotonicities in sev¬ 
eral thermodynamic quantities, including Szz and den¬ 
sity, across free interfaces of the films that do not un¬ 
dergo surface freezing. We do not detect any such non¬ 
monotonicity in relaxation time profiles. Yet, we are un¬ 
able to rule out its possibility for the following reasons. 
First of all, the uncertainties associated with computing 
dynamical quantities such as relaxation times are gener¬ 
ally larger than the uncertainties in estimating thermo¬ 
dynamic quantities such as density. After all, the former 
are computed from auto-correlation functions while the 
later are obtained from simple time averaging. Secondly, 
the convention used for defining ^-dependent autocorre¬ 
lation functions will inevitably lead to some mixing be¬ 
tween neighboring slices as the molecules that contribute 
to the autocorrelation function of a particular 2 ; slice can 
leave and re-enter that slice. Such a procedure can mask 
the existence of subtle spatial differences in relaxation 
time, a situation that is completely plausible in the pre¬ 
frozen interfaces considered in this work. 

The dynamical features of a film in the vicinity of 
a substrate depend on e^. Like the free interface, 
the surface-frozen monolayers close to loose substrates 
demonstrate a deceleration of dynamics and a decoupling 
of translational and rotational degrees of freedom. This 
slow-down happens in the vicinity of a sticky substrate 
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FIG. 10: Profiles of normal and lateral stress for different temperatures and es values. 
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FIG. 11: Gorrelation between the profiles of atomic density and normal stress for a film at (a) T = 260 K, = 0.5, and (b) 
T = 240 K, es = 1.0. The region shaded in light blue has long-range translational order. 


as well. However, it is not clear whether translational 
and rotational degrees of freedom decouple as the cor¬ 
responding relaxation times cannot be computed in the 
timescale of our simulations. For the films that do not 
undergo surface freezing, however, two distinct dynami¬ 
cal regimes are observed. For sticky walls {es = 3.0), the 
dynamics is decelerated close to the substrate, even in 
amorphous regions of the film. Indeed, relaxation times 
can be as much as two orders of magnitude larger than 
the corresponding bulk values. In the vicinity of loose 
substrates {es = 0.5), however, the dynamics is acceler¬ 
ated, as in atomic thin films. For es = 1.0, the relaxation 
time profiles are more oscillatory with the overall dynam¬ 


ics becoming only slightly faster with respect to the bulk. 
This is consistent with the behavior observed in atomic 
thin films that also demonstrate two distinct dynamical 
regimes in the vicinity of sticky and loose substrates. 

Fig. 14 depicts the temperature dependence of bulk 
relaxation times for es = 0.5. Those are obtained from 
averaging the relaxation time profiles of Figs. 12 and 13 
over 2 :[nm] G [2.5, 7.5]. An identical behavior is observed 
for other es values, with the results not shown for con¬ 
ciseness. All relaxation times demonstrate an Arrhneius- 
type dependence on temperature, with no sign of fragility. 
This is not unusual considering the earlier viscosity mea¬ 
surements for n-octane [86, 87], which predict a glass 
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FIG. 12: Profiles of atomic and molecular translational relaxation time for different temperatures and es values. 
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FIG. 13: Profiles of bond and molecular rotational relaxation time for different temperatures and values. 





FIG. 14: Temperature dependence of bulk relaxation times 
for es = 0.5. 


transition temperature of = 85 K, markedly lower 
than the temperatures considered here. Note that atomic 
translational relaxation times are always lower than their 
molecular counterparts, due to the relative ease of re¬ 
laxing atomic degrees of freedom. It is, however, not 
possible to systematically compare the bond and molec¬ 
ular orientational relaxation times, since they correspond 
to relaxation features of different Legendre polynomials. 
Therefore, the fact the bond relaxation times are higher 
does not have any physical meaning. Similarly, it is not 
meaningful to compare translational and orientational re¬ 
laxation times, as they have also been computed from 
decays of different autocorrelation functions. 

In our earlier work [5], we established a close corre¬ 
lation between lateral relaxation time and lateral stress 
profiles in amorphous regions of atomic thin films. This is 
consistent with what happens in simple liquids in which 
diffusivity increases upon decreasing pressure [88]. We 
find a similar correlation between the profiles of atomic 
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FIG. 15: Correlation between the profiles of atomic translational relaxation time and lateral stress for a film at (a) T = 260 K, 
es = 0.5, and (b) T = 240 K, es = 1.0. The region shaded in light blue has long-range translational order. 


translational relaxation time and lateral stress in octane 
films. Fig. 15 shows such a correlation for two represen¬ 
tative films, with the peaks and valleys of Ttr,a and lateral 
stress closely following one another in the amorphous re¬ 
gion (determined from RDF). It is interesting to note 
that such a correlation can also be rationalized by the 
fact that the compressibility of a simple fluid increases 
upon decreasing pressure, and this enhances structural 
relaxation by facilitating local density fluctuations. 

IV. CONCLUSIONS 

In this work, we study the effect of substrate on ther¬ 
modynamic and kinetic anisotropies in n-octane thin 
films. We observe complete freezing at temperatures 
below 232.5 ± 2.5 K. For 235 K < T < 250 K, how¬ 
ever, a frozen monolayer emerges at vapor-liquid inter¬ 
faces, as well as in the vicinity of loosely attractive 
substrates. Such frozen monolayers correspond to large 
peaks in several thermodynamic and kinetic properties, 
such as density, potential energy, nematic order param¬ 
eter, and translational and rotational relaxation times. 
Also, translational and rotational degrees of freedom are 
decoupled in the monolayer due to long-range orienta¬ 
tional ordering of the n-octane molecules. 

At higher temperatures, interfacial freezing only occurs 
in the vicinity of sticky substrates. At vapor-liquid inter¬ 
faces, and close to loose substrates, only a weak propen¬ 
sity is observed for the molecules to ’pre-freeze‘ and align 
perpendicular to the interface. Such a propensity mani¬ 
fest itself in mild peaks in atomic and molecular density 
and f{0,z) and Szz orientational order parameters. 

The amorphous regions of the films are stratified in the 
vicinity of substrates, with the extent of stratification in¬ 
creasing upon decreasing temperature, or increasing the 


interaction parameter e^. Density oscillations are min¬ 
imal close to loose substrates, and the interface almost 
resembles a vapor-liquid interface. Similar to density, os¬ 
cillations are observed in other thermodynamic quantities 
such as potential energy, dive (as a result of energy mini¬ 
mization), stress, and orientational order parameters. In 
the vicinity of a substrate, we confirm the existence of 
correlation between atomic density and normal stress in 
amorphous regions of the film, consistent with our ob¬ 
servations in atomic thin films. We also observe two dis¬ 
tinct dynamical regimes at the solid-liquid interface. In 
the vicinity of loose substrates, dynamics is accelerated, 
while a sticky substrate decelerates dynamics in its vicin¬ 
ity. Finally, we are able to establish a correlation between 
lateral atomic translational relaxation times and lateral 
stress in amorphous regions of the films. 

Short-chain n-alkanes are known to be poor glass form¬ 
ers due to their strong propensity to crystallize [89] . Our 
findings suggest that the vapor deposition process that 
yields ultrastable glasses for many materials is unlikely 
to be very successful in the case of n-alkanes. This is 
not only due to the propensity of alkane films to undergo 
surface freezing, but also because of the fact that the free 
interface does not create a more accessible potential en¬ 
ergy landscape even when the alkane film is amorphous. 
It has, indeed, been recently demonstrated that there are 
materials that cannot form a ultrastable glass upon va¬ 
por deposition [90], unlike earlier suggestions that this 
process might be universal [91]. This has led to specula¬ 
tions that a correlation might exist between the fragility 
of a liquid and its ability to form an ultrastable glass [92] . 
This correlation is not clear-cut as some relatively strong 
liquids still form ultrastable glasses upon vapor deposi¬ 
tion [93]. There might, however, be other dimensions to 
this problem, namely the presence of a mobile interfa¬ 
cial region, as well the accessibility of deeper minima of 
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the potential energy landscape at the interface. If this 
picture is accurate, n-octane must not yield ultrastable 
glasses upon physical vapor deposition. It is therefore a 
worthwhile experimental endeavor to investigate vapor- 
deposited glasses of n-alkanes, and other molecules with 
aliphatic chains to determine whether this picture is ac¬ 
curate. 

It is necessary to emphasize that, according to this pic¬ 
ture, elevated orientational ordering at a free interface 
will not necessarily make a material a poor ultrastable 
glass former. Indeed, experimental [30-32] and compu¬ 
tational [50] studies of ultrastable glasses have confirmed 
the existence of such ordering at vapor-liquid interfaces 
of good ultrastable glass formers. 

It is necessary to emphasize that the microstructure 
of the liquid in the vicinity of the solid-liquid interface 
might depend on the structure of the substrate. For in¬ 
stance, the orientational order that would emerge close 
to a sticky structureless substrate will be different from 
what is observed here, as reported in Ref. [94] . Since the 
liquid that is close to a low-e^ substrate is structurally 
very similar to the liquid at the free interface, we expect 
that such differences will only be important at large val¬ 
ues of e^. Even then, the observed deceleration of the 
dynamics (with respect to bulk) is not expected to dis¬ 
appear if a structureless sticky substrate is utilized. It is, 
however, interesting to study the combined effect of 
and substrate corrugation on the microstructure of the 


frozen monolayers, as well as the thermodynamic and ki¬ 
netic anisotropies in the film. 

There are other interesting questions to be asked about 
alkane films besides their ability or lack thereof to form 
ultrastable glasses. One such question is the role of sub¬ 
strate in inducing- or suppressing- surface freezing. This 
can be systematically addressed by performing rigorous 
free energy and/or rate calculations. Finally, the ultra- 
thin films of long-chain alkanes can be studied to in¬ 
vestigate the possible interplay between solid-liquid and 
vapor-liquid interfaces. All these topics can be the sub¬ 
ject of further studies. 
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